; docformat = 'rst'
;
;+
;
; :Purpose:
;   Construct growth-constraint matrices
;   
; :Inputs:
;
; :Outputs:
;   Structure containing:
;     'GROWTH': expectation value of growth constraint
;     'D': Time-based differencing operator
;     'DTS1': transpose(D)##S^(-1)
;     'DTS1D': transpose(D)##S^(-1)##D
;     'WH': location of time-based differences
;     'WH2D': location of time-based differences in 2D (for square matrices)
;     
;     'D_L': Spatial differencing operator
;     'DTS1_L': transpose(D_L)##S_L^(-1)
;     'DTS1D_L': transpose(D_L)##S_L^(-1)##D_L
;     'WH_L': location of latitude-based differences
;     'WH2D_L': location of latitude-based differences in 2D (for square matrices)
;
; :Example::
;
; :History:
; 	Written by: Matt Rigby, MIT, Jun 20, 2013
;
;-
pro a2i_growth, parameters, x, emissions, measurements

  x_growth=!null
  x_growth_error=!null
  x_growth_wh=!null
  
  x_growth_L=!null
  x_growth_error_L=!null
  x_growth_wh_L=!null

  for pi=0, parameters['NPOLLUTANTS']-1 do begin

    nYears_species=parameters['NYEARS'] - measurements['START', pi]/12L

    ;Initial conditions spatial constrain
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    if parameters['IC_GROWTH', pi] then begin
      for bi=0, 1 do begin
        wh_start=where(x['POLLUTANT'] eq pi and x['ID'] eq 'I' and x['BOX'] ge bi*4 and x['BOX'] lt (bi+1)*4, count_start)
        wh_end=where(x['POLLUTANT'] eq pi and x['ID'] eq 'I' and x['BOX'] ge (bi+1)*4 and x['BOX'] lt (bi+2)*4, count_end)
        if count_start gt 0 then begin
          x_growth_wh_L=[x_growth_wh_L, [[wh_start], [wh_end]]]
          x_growth_L=[x_growth_L, replicate(0., count_start)]
          x_growth_error_L=[x_growth_error_L, replicate(parameters['ERROR_IC', pi], count_start)]
        endif
      endfor
      ;Also add latitudinal cross-section at surface
      wh_start=where(x['POLLUTANT'] eq pi and x['ID'] eq 'I' and x['BOX'] ge 0 and x['BOX'] lt 3, count_start)
      wh_end=where(x['POLLUTANT'] eq pi and x['ID'] eq 'I' and x['BOX'] ge 1 and x['BOX'] lt 4, count_end)
      x_growth_wh_L=[x_growth_wh_L, [[wh_start], [wh_end]]]
      x_growth_L=[x_growth_L, replicate(0., count_start)]
      x_growth_error_L=[x_growth_error_L, replicate(parameters['ERROR_IC', pi], count_start)]
    endif

    ;Emissions growth constraint
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    if (parameters['EMISSIONS_ESTIMATE', pi]) and $
      (parameters['EMISSIONS_GROWTH', pi]) and $
      (nYears_species gt 1) then begin
      
      ;Emissions time constraint
      for bi=0, 3 do begin
        wh=where(x['POLLUTANT'] eq pi and x['ID'] eq 'Q' and x['BOX'] eq bi, count)
        if count gt 0 then begin
          wh_start=wh[0:-1 - parameters['EMISSIONS_FREQUENCY', pi]]
          wh_end=intarr(n_elements(wh_start))
          for n=0, n_elements(wh_start)-1 do begin
            whYear=where(x['TI', wh] eq x['TI', wh[n]]+12, countYear)
            if countYear ne 1 then stop
            wh_end[n]=wh[whYear]
          endfor
          x_growth_wh=[x_growth_wh, [[wh_start], [wh_end]]]
          error=reform(emissions['EMISSIONS_GROWTH_ERROR', pi, bi, *])
          error=interpol(error, $
            findgen(n_elements(error))/12., $
            findgen(count - parameters['EMISSIONS_FREQUENCY', pi])/float(parameters['EMISSIONS_FREQUENCY', pi]) + 0.5)
        endif else begin
          message, 'SOMETHING WEIRD HERE'
        endelse
        
        if n_elements(error) ne n_elements(wh_end) then message, 'SOMETHING WRONG WITH GROWTH CALCULATION'
        x_growth_error=[x_growth_error, error]
        x_growth=[x_growth, replicate(0., n_elements(wh_end))]

      endfor  ;box

      ;Emissions spatial constraint
      for bi=0, 2 do begin
        wh_start=where(x['POLLUTANT'] eq pi and x['ID'] eq 'Q' and x['BOX'] eq bi, count_start)
        wh_end=where(x['POLLUTANT'] eq pi and x['ID'] eq 'Q' and x['BOX'] eq bi+1, count_end)
        if count_start gt 0 then begin
          x_growth_L=[x_growth_L, replicate(0., count_start)]
          case parameters['FUNCTION_Q', pi] of
            'X': begin
              x_growth_error_L=[x_growth_error_L, $
                replicate(parameters['EMISSIONS_LATITUDINAL_GROWTH_ERROR', pi], count_start)]
            end
            'LN(X)': begin
              x_growth_error_L=[x_growth_error_L, $
                replicate((parameters['EMISSIONS_LATITUDINAL_GROWTH_ERROR', pi]), count_start)]
            end
          endcase
          x_growth_wh_L=[x_growth_wh_L, [[wh_start], [wh_end]]]
        endif
      endfor  ;box
      
    endif

  endfor  ;pollutant

  ;Model parameters
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  for paramI=0, 1 do begin
    case paramI of
      0: begin
       param='PT'
       nBox=17
       uncertainty=replicate(parameters['ERROR_T'], nBox)
       do_parameter=parameters['T_GROWTH']
      end
      1: begin
        param='PO'
        nBox=4
        uncertainty=parameters['ERROR_OH']
        do_parameter=parameters['OH_GROWTH']
      end
    endcase
    
    if do_parameter then begin
      for bi=0, nBox-1 do begin
        wh=where(x['ID'] eq param and x['BOX'] eq bi, count)
        if count gt 0 then begin
          wh_start=wh[0:-5] ;Both parameters are seasonal
          wh_end=intarr(n_elements(wh_start))
          for n=0, n_elements(wh_start)-1 do begin
            whYear=where(x['TI', wh] eq x['TI', wh[n]]+12, countYear)
            if countYear ne 1 then stop
            wh_end[n]=wh[whYear]
          endfor
          x_growth=[x_growth, replicate(0., n_elements(wh_end))]
          x_growth_error=[x_growth_error, replicate(uncertainty[bi], n_elements(wh_end))]
          x_growth_wh=[x_growth_wh, [[wh_start], [wh_end]]]
        endif
      endfor  ;box
    endif
  endfor

  ;Construct matrices
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  for gi=0, 1 do begin
    
    case gi of
      0: begin
        wh=x_growth_wh
        error=x_growth_error
        growth=x_growth
        latStr=''
      end
      1: begin
        wh=x_growth_wh_L
        error=x_growth_error_L
        growth=x_growth_L
        latStr='_L'
      end
    endcase
    
    if n_elements(wh) gt 0 then begin
      D_temp=fltarr(x['NSTATE'], x['NSTATE'])
      D_temp[wh[*, 0], wh[*, 0]]=-1.
      D_temp[wh[*, 1], wh[*, 0]]=1.

      S=replicate(0., x['NSTATE'])
      S[wh[*, 0]]=error^2
      S1=fltarr(x['NSTATE'])
      whS=where(S gt 0.)
      S1[whS]=1./S[whS]

      x_growth_error_temp=replicate(0., x['NSTATE'])
      x_growth_error_temp[wh[*, 0]]=error
  
      x_Growth_temp=replicate(0., x['NSTATE'])
      x_Growth_temp[wh[*, 0]]=growth
      
      ;DTS1=transpose(D)##diag_matrix(S1)
      DTS1_temp=transpose(D_temp)*rebin(S1, [x['NSTATE'], x['NSTATE']])
      DTS1D_temp=DTS1_temp##D_temp
      
      ;Extract only growth-related state-vector elements
      xWhG_temp=fltarr(x['NSTATE'])
      xwhG_temp[wh]=1
      xWhG_temp=where(xWhG_temp, nG)

      if nG gt 0 then begin
      
        xWhG2D_temp=lonarr(nG^2)
        for n=0L, nG-1 do xWhG2D_temp[n*nG:(n+1)*nG-1]=xWhG_temp[n]*x['NSTATE'] + xWhG_temp
        
        x_growth_temp=x_Growth_temp[xWhG_temp]
        x_growth_error_temp=x_Growth_error_temp[xWhG_temp]
        
        D_temp=D_temp[xwhG_temp, *]
        D_temp=D_temp[*, xwhG_temp]
        
        DtS1_temp=DtS1_temp[xwhG_temp, *]
        DTS1_temp=DTS1_temp[*, xwhG_temp]
        
        DtS1D_temp=DtS1D_temp[xwhG_temp, *]
        DTS1D_temp=DTS1D_temp[*, xwhG_temp]
        
      endif else begin
  
        xWhG_temp=-1
        xWh2D_temp=-1
        x_growth_temp=-9999
        D_temp=-9999
        DTS1_temp=-9999
        DTS1D_temp=-9999
        
      endelse

      x['GROWTH' + latStr]=x_growth_temp
      x['GROWTH_ERROR' + latStr]=x_growth_error_temp
      x['D' + latStr]=D_temp
      x['DTS1' + latStr]=DTS1_temp
      x['DTS1D' + latStr]=DTS1D_temp
      x['GROWTH_WH' + latStr]=xWHG_temp
      x['GROWTH_WH2D' + latStr]=xWHG2D_temp
      
    endif
    
  endfor


;  for bi=0, 2 do begin
;    wh1=where(x['ID'] eq 'Q' and x['BOX'] eq bi, count1)
;    wh2=where(x['ID'] eq 'Q' and x['BOX'] eq bi+1, count2)
;    if count1 gt 0 then begin
;      for n=0, count1-1 do begin
;        D_L[wh1[n], wh1[n]]-=1.
;        D_L[wh2[n], wh1[n]]=1.
;      endfor
;    endif
;    S_L[wh1]=10.^(2.)
;    growth_wh_latitude=[growth_wh_latitude, [[wh1], [wh2]]]
;  endfor

;  ;latitudinal growth constraint
;  D_L=fltarr(x['NSTATE'], x['NSTATE'])
;  S_L=replicate(0., x['NSTATE'])  
;  S_L1=fltarr(x['NSTATE'])
;  wh=where(S_L gt 0.)
;  S_L1[wh]=1./S_L[wh]
;  ;DTS1_L=transpose(D_L)##diag_matrix(SL1)
;  DTS1_L=transpose(D_L)*rebin(S_L1, [x['NSTATE'], x['NSTATE']])
;  DTS1D_L=DTS1_L##D_L
;
;  
;  ;Extract only growth-related state-vector elements (Latitudinal growth constraint)
;  xWhG_L=fltarr(x['NSTATE'])
;  xwhG_L[growth_wh_latitude]=1
;  xWhG_L=where(xWhG_L, nG_L)
;  
;  xWhG2D_L=lonarr(nG_L^2)
;  for n=0L, nG_L-1 do xWhG2D_L[n*nG_L:(n+1)*nG_L-1]=xWhG_L[n]*x['NSTATE'] + xWhG_L
;  
;  D_L=D_L[xWhG_L, *]
;  D_L=D_L[*, xWhG_L]
;  
;  DtS1_L=DtS1_L[xWhG_L, *]
;  DTS1_L=DTS1_L[*, xWhG_L]
;  
;  DtS1D_L=DtS1D_L[xWhG_L, *]
;  DTS1D_L=DTS1D_L[*, xWhG_L]

;  return, {GROWTH:xGrowth, D:D, DTS1:DTS1, DTS1D:DTS1D, WH:xWhG, WH2D:xWhG2D, $
;    D_L:D_L, DTS1_L:DTS1_L, DTS1D_L:DTS1D_L, WH_L:xWhG_L, WH2D_L:xWhG2D_L}

end